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We introduce a cosmological model based on the normal branch of DGP braneworld gravity with 
a smooth dark energy component on the brane. The expansion history in this model is identical to 
ACDM, thus evading all geometric constraints on the DGP cross-over scale r c . This well-defined 
0^ , model can serve as a first approximation to more general braneworld models whose cosmological 

■ solutions have not been obtained yet. We study the formation of large scale structure in this model in 
the linear and non- linear regime using N-body simulations for different values of r c . The simulations 
use the code presented in [1] and solve the full non-linear equation for the brane-bending mode in 

. conjunction with the usual gravitational dynamics. The brane-bending mode is attractive rather 

than repulsive in the DGP normal branch, hence the sign of the modified gravity effects is reversed 
compared to those presented in [1]. 

We compare the simulation results with those of ordinary ACDM simulations run using the 
same code and initial conditions. We find that the matter power spectrum in this model shows a 
characteristic enhancement peaking at k ~ 0.7/i/Mpc. We also find that the abundance of massive 
halos is significantly enhanced. Other results presented here include the density profiles of dark 
matter halos, and signatures of the brane-bending mode self-interactions (Vainshtein mechanism) 
in the simulations. Independently of the expansion history, these results can be used to place 
, constraints on the DGP model and future generalizations through their effects on the growth of 

T ■ cosmological structure. 
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■ I. INTRODUCTION 

C/2 . 

Braneworld scenarios with infinite extra dimensions have received much interest in recent years, as a possible 
, explanation for the observed accelerated expansion of the universe [2, 3], as well as a way to tackle the cosmological 
J> ■ constant problem [4, 5]. The simplest such model is the Dvali-Gabadadze-Porrati (DGP) model [6]. In this model, 
| matter lives on a four-dimensional brane embedded in five-dimensional Minkowski space. The gravity action consists 
C*~) . of a five-dimensional Einstein-Hilbert term plus an ordinary, four-dimensional term localized on the brane. Gravity 
thus is five-dimensional on large scales, and reduces to four-dimensional General Relativity on small scales. The 
, transition is given by the cross-over scale r c , which is defined by the ratio of five- and four-dimensional gravitational 
constants: r c = G ( ^/2G^\ 

There are two branches of homogeneous and isotropic solutions in this model, corresponding to the two possible 
ways of embedding the brane in (asymptotic) five-dimensional Minkowski space. The self-accelerated branch (sDGP) 
has received attention since it leads to a late-time accelerated expansion of the universe without any dark energy 
or cosmological constant, if r c is of order the present Hubble horizon c/Hq. However, this branch of the model is 
plagued by a ghost instability when perturbed around the de Sitter solution [7-9] . In addition, the expansion history 
predicted by the self-accelerating DGP model [2] does not appear to fit observations (e.g., [10]). The other, so-called 
normal branch of the theory (nDGP), does not have a ghost instability. However, it does not lead to self-acceleration, 
so that it is necessary to add a cosmological constant (tension) on the brane [11-14], or, more generally a form of 
stress energy with negative pressure. 

Hence, we are interested in generalizing the normal-branch DGP model in a way that allows it to pass expansion 
history constraints. Much work is ongoing to extend the DGP braneworld scenario [5, 15, 16], and the proposed models 
are generally expected to exhibit an expansion history close to ACDM, while the phenomenology of the modification 
to GR is expected to be similar to the normal branch DGP model. However, full cosmological solutions have not been 
obtained yet. Braneworld-inspired parametrized generalizations of DGP have been proposed in [17, 18]. On the one 
hand, this approach might well capture features of a yet-unknown full generalized braneworld model. On the other 
hand, these parametrizations do not constitute a well-defined theory. 
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Here, we instead opt for considering a well-defined, albeit somewhat contrived model based on DGP, which nev- 
ertheless exhibits some of the features of a generalized braneworld model (see Section II D for a discussion). Our 
model consists of normal-branch DGP gravity, with a general, smooth quintessence-type dark energy component on 
the brane tuned so that the resulting expansion history is identical to a ACDM model. This model which we call 
nDGP+DE thus satisfies all geometric constraints from, e.g. the cosmic microwave background (CMB) and Super- 
novae, and leaves the cross-over scale r c as an almost free parameter of the model. Furthermore, the r c — > oo limit of 
the theory is precisely General Relativity with a cosmological constant. Our main aim in this paper is to study the 
growth of structure in nDGP+DE in the linear and non-linear regime. 

The evolution of linear perturbations in DGP has been studied in [19-22]. On scales smaller than the horizon 
and the cross-over scale, the model can be described as an effective scalar-tensor theory. The massless scalar degree 
of freedom corresponds to displacements of the brane and is called the brane-bending mode. In linear perturbation 
theory, the effect of the brane-bending mode is simply to rescale the effective gravitational constant for the dynamical 
potential (the time-time piece of the metric), while the geodesies of photons (i.e. gravitational lensing) are not affected. 
In the self-accelerating branch of DGP, the brane-bending mode is repulsive, weakening the effective gravitational 
force, while it is attractive in the normal branch. 

As all viable modified gravity models, DGP contains a non-linear mechanism to restore General Relativity in 
high-density environments. This is achieved by a strong self-coupling of the brane-bending mode, which becomes 
effective within a characteristic scale, the so-called Vainshtein radius [23-25]. The Vainshtein radius depends on the 
mass considered as well as its configuration. While the prefactor of the self-coupling is model-dependent, the form 
of the coupling is generic and is expected to be universal to braneworld models [15, 22]. For typical structures in 
the Universe, the Vainshtein radius is of cosmological scale. Hence it is necessary to follow the brane-bending mode 
and its self- interactions when studying the formation of structure in the Universe. In [1], we presented simulations 
of the self-accelerating DGP model which self-consistently solve for the brane-bending mode (see also [26]). Here, we 
present simulations of two nDGP+DE models, extending the studies of non-linear structure formation presented in 
[1] to the normal branch of DGP, where the effects of modified gravity have the opposite sign. 

N-body simulations of normal-branch braneworld models have previously been presented in [18]. The key difference 
to our simulations is that while we solve the full brane-bending mode equation, [18] employed an approximation where 
the modified forces are parametrized by an effective density-dependent gravitational constant G c s{8p). Unfortunately, 
as already pointed out in [18], this approximation does not recover the correct large-distance behavior for the brane- 
bending mode, so that in principle artefacts of the approximation might appear even on large scales. Furthermore, the 
approximation will become worse the higher the resolution of the simulation is. We present a quantitative comparison 
of our simulations with the G e e(5) approximation in Appendix A. 

In addition to the matter power spectrum and halo mass function, we show results on the density profiles of halos, 
and signatures of the Vainshtein mechanism in the simulation results. The results presented in this paper together 
with [1] can serve as a basis for modeling non-linear structure formation in DGP. A model of some of the results 
presented here and in [1] in the context of spherical collapse and the halo model will be the subject of a forthcoming 
paper. Such a model can then be used to constrain DGP and more generalized braneworld models using the wealth 
of large scale structure observations available (see [27] for a first attempt in the case of f(R) gravity). 

The paper is structured as follows. In Section II, we describe the nDGP+DE models in detail, including the 
evolution of linear perturbations and the Vainshtein mechanism. We describe the simulations in Section III, and 
present the results in Section IV. We conclude in Section V. 

II. DGP MODELS 

As a gravitational framework for our nDGP+DE cosmology, we choose the normal branch of the DGP model [2]. 
However, we add a smooth, quintessence-type dark energy on the brane whose equation of state is adjusted precisely 
to cancel the unwanted effects of the extra dimension on the expansion history, yielding an expansion history identical 
to ACDM. Hence, independently of the expansion history, the cross-over scale r c is a free parameter which is not 
constrained anymore to be of order Hq 1 in this model. The r c — > oo limit of this theory is General Relativity with a 
cosmological constant. 

We will consider flat cosmologies throughout, and consider two cosmologies of the nDGP+DE type: one with 
r c = 500 Mpc ("nDGP-1"), and one with r c = 3000 Mpc ("nDGP-2"). We will also compare our results with the 
cosmological simulations of the self-accelerating DGP model without dark energy ("sDGP") presented in [1]. 
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FIG. 1: Behavior of the dark energy component in the nDGP+DE models: density pde and p m , as well as the resulting effective 
PA {left panel); equation of state hide of the dark energy (right panel). The model parameters are defined in Tab. I. 



A. Background evolution 



We start with the modified Friedmann equation for the normal branch of DGP (e.g. [11, 14]): 

_ H(a) _ 
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where Hq is the Hubble rate a/ a today, and 
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Here p m and pr>R stand for the background energy density of matter and smooth dark energy, respectively. We can 
now simply equate Eq. (2.1) with the Friedmann equation in GR for a flat ACDM universe (with il\ = 1 — Q m ): 



E 



ACDM 



(a) = \/Vl m a 3 + £!a- 



(2.3) 



and solve for the dark energy density. Note that in this model, one would infer the correct value of Q m from geometric 
and early-Universe tests, such as BBN, CMB, SN, and Hq measurements, which rely on the assumption of the 
Friedmann equation Eq. (2.3) 1 . This greatly simplifies the interpretation of our results below. We then obtain the 
following dark energy density: 



Pde (a) = Pa,o 



1 - 1 



Note that this energy density is positive definite. In particular, we can derive the following two limits: 



Pde (a) = Pcr,o 



Q\ + 2f2 rc 



(2.4) 



(2.5) 



Of course, the interpretation of A as a cosmological constant on the other hand would be very misguided. 
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FIG. 2: Left panel: The f3 function Eq. (2.9) in the nDGP+DE models. We also show the function for the self-accelerating 
DGP model considered in [1] (note that /3 < for that model). Right panel: Linear growth factor D(a)/a as a function of a, 
where D(a) is normalized to a in the matter-dominated regime, for ACDM and nDGP+DE models in the quasi-static regime 
(k > 0.01/i/ M P c )- 



This shows that the equation of state u>de = Pde/pde is always greater than -1: it approaches — 1/2 for s<l, and 
— 1 for a > 1 (Fig. 1). For r c — ► oo, i.e. f2 rc — > 0, pde = ^APcr.o =const., confirming the ACDM limit. Note that in 
the matter-dominated regime, the dark energy density increases with decreasing cross-over scale r c (i.e., increasing 
fi rc )j while at late times, when the DGP term dominates, it becomes independent of r c and only depends on the 
effective ACDM parameters. Fig. 1 shows the evolution of the dark energy density and the equation of state for the 
two models we have simulated, nDGP-1 (r c = 500 Mpc) and nDGP-2 (r c = 3000 Mpc). 



B. Linear cosmological perturbations 

We now consider the evolution of cosmological perturbations in our nDGP+DE model. We work in conformal 
Newtonian gauge and write the metric as: 

ds 2 = -[1 + 2^{x,t)]dt 2 + [1 + 2$(x,i)]a 2 (i)dx 2 (2.6) 

Since we assume that the additional dark energy component is smooth, the growth of structure in this model proceeds 
in the same way as in the normal branch with brane tension [12-14], except that the expansion history is altered. 
In particular, on small scales k > 0.01 /i/Mpc, time derivatives may be neglected with respect to spatial derivatives, 
which we call the quasi-static regime. If, in addition, k ^S> r" 1 , DGP reduces to an effective scalar-tensor theory, 
where the brane-bending mode ip plays the role of the additional degree of freedom. The linearized brane-bending 
mode equations then lead to a scale-independent but time-dependent modification of the linear growth equation for 
matter perturbations, 5 = (p — p m )/p m (e.g. [20, 25, 28]): 

4r* = i*GU<*)T'„A (2-7) 
<W«)=g(i+3^)), (2-8) 
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where dots denote time derivatives, G is Newton's constant, and, in the normal branch of DGP: 

( 3(a) = l + 2H(a)r c (^ + J^ ) j. (2.9) 

The function /3 thus quantifies the departure from GR, which is restored when (3 3> 1 at early times (Fig. 2, left 
panel). In the normal branch, (3 > 1 and hence graviational forces are enhanced by up to 1/3. Conversely, in the 
self-accelerating branch /3 < — 1 and forces are suppressed. Since we only add an additional smooth stress-energy 
component and do not modify the gravity sector, the changes to standard DGP only come in through the expansion 
history in Eq. (2.9). 

Fig. 2 (right panel) shows the matter growth factor D(a) = 5(a) /8{ai) in the quasi-static regime as a function 
of a. The normalization is such that D(a) = a in the matter-dominated regime. As in GR+DE, the growth factor 
is independent of scale in the quasi-static regime (k 3> r~ ). In this case, the growth is always enhanced in the 
nDGP+DE models with respect to ACDM. On such scales much smaller than r c , the five-dimensional nature of 
gravity is not important, and the effect of the attractive brane-bending mode is the cause for the differences to 
ACDM. 

More generally, the effect of the brane-bending mode ip on the metric potentials 'J, $ on sub- horizon and sub-r c 
scales is given by: 

* = + (2.10) 
$ = -*jv + ^- (2.11) 

Here, we have defined the "Newtonian" potential via the usual Poisson equation, 

k 2 

— V N - 47rGp m <5. (2.12) 

Note that the propagation of light is not directly affected in DGP, as is usually the case in simple scalar-tensor theories: 
in the quasi-static regime, the lensing potential $_ = (^ — ( f>)/2 reduces to <3>Ar. In other words, the relation between 
<&_ determining the geodesies of photons and the matter overdensities is unchanged from GR [28]. 

Finally, on larger scales approaching the horizon, time derivatives as well as the effects of the extra dimension 
become important. Here, one can use the parametrized post-Friedmann approach [29], calibrated on calculations 
using the dynamical scaling solution [20], or the full numerical solution [21]. This approach was also used in [14] to 
place constraints on DGP models with cosmological constant (brane tension), and we adopt the parameters for the 
normal branch given in [14]. Since we are considering models with r c smaller than the present horizon scale, care has 
to be taken on scales significantly larger than r c where the PPF parametrization has not been calibrated with the full 
5D calculation. Such large scales are never relevant for the simulations presented here, however. 

Like all N-body simulations, our simulations work in the quasi-static approximation. The full DGP equations 
contain non-local terms involving Vk 2 /r c which could be relevant for small r c . We verify that neglecting these terms 
has a negligible impact on the scales probed by the simulations in Section III. 



C. Vainshtein mechanism 



If the modified force law given by Gu n ^ G was valid on all scales, the DGP model would predict order unity 
deviations from the GR values of post-Newtonian parameters in the Solar System. However, the brane-bending mode 
<p has self-interactions which suppress its value in high-density regions. In the quasi-static regime, the equation of 
motion can be written as: 

VV+ 3^ [(VV)2 _ ( V * V ^)( V * V -V)] = ^0^8p, (2.13) 

where 5p = ~p m 5. The motion of particles is then governed by the dynamical potential as in Eq. (2.10), which 
receives a contribution from <p. Linearizing Eq. (2.13) together with Eq. (2.10) then yields Eqs. (2.7)-(2.8). For 
spherically symmetric masses, one can derive an analytic solution to the field equation. Assuming a mass of constant 
density with radius R (top hat), one obtains the following gravitational acceleration: 

5 = SN + ~=S N (r)[l + A(r)}, (2.14) 
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FIG. 3: Overdensity threshold for the tp self-interactions, defined as inverse of the non-linearity parameter g [Eq. (2.17)], as a 
function of scale factor in the different DGP cosmologies. 



A(r) 



2 

3/3 



v R , 



1 



r > R 



r < R. 



(2.15) 



where g N is the Newtonian acceleration of the spherical mass, and r* denotes a characteristic scale of the solution, 
the Vainsktein radius, determined by r c and the Schwarzschild radius of the mass r s : 



9/3 2 



(2.16) 



Close to the mass (r <C r»), force modifications are suppressed by ~ (r/r») 3 / 2 . At very large distances, r ^> r», A(r) 
approaches the constant 1/(3/3), which exactly matches the linear solution, Eq. (2.8). Note that for the Sun and 
r c = 3000 Mpc, r„ ~ 75 pc, while for the Milky way r* ~ Mpc. Hence, for r c ~ Gpc the non-linear interactions in 
Eq. (2.13) are expected to be important on cosmological scales. Thus it is crucial to self-consistently solve for the 
non-linearities in tp in conjunction with the evolution of structure in the Universe. 

Another limiting case in which Eq. (2.13) is easily solvable is a plane wave density field: in this case, the two 
non-linear terms cancel, and one recovers the linearized solution. One might expect cosmological structure, with its 
sheets and filaments as well as virialized structures, to lie somewhere in between these limiting cases. Hence, the 
precise structure of the non-linearity in Eq. (2.13) is important for structure formation. 

An order-of- magnitude estimate of where the coupling of the brane-bending mode becomes strong can be obtained 
as follows. We take the ratio of the linear and non-linear ip terms in Eq. (2.13) and evaluate them for the linear 
solution tp l : 



2^2 



3/3 a- 



■(vV) 2 / vV 



9/3 2 9 3/3 2 



■a 



■5 = g5, 



(2.17) 



One can then estimate that wherever S > 5th = the non-linear interactions of tp are important (see Fig. 3). The 
quantity g defined here is the non-linearity parameter introduced in Sec. IV C of [22]. At early times, g — * 1/3. By 
z = 0, it decreases to 10 -3 and 10 -2 in the nDGP-l and nDGP-2 cosmologies, respectively, and ~ 0.1 in sDGP. 
Today then, the non-linearity criterion will be satisfied for S > 0(1000) in the case of nDGP-l, and O(100) in the 
case of nDGP-2, which will only be the case within the inner regions of dark matter halos. In the case of sDGP, 
<5th w 10 today, which makes the non-linearities important in a significant fraction of the Universe [1]. 
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D. Comparison with parametrized braneworld models 



The nDGP+DE model presented here is a full consistent (albeit somewhat contrived) theory. It is worth comparing 
this model with the parametrized braneworld-inspired models introduced in [17, 18]. On sub-horizon scales, k 3> 
aH, r c , these effective models yield the same equations as nDGP+DE (Section II B and II C), except for a modified 
function (3: 



For a = 1/2, this reduces to Eq. (2.9). Hence, the a = 1/2 model of [18] is identical to our model on small scales. 
Although the a = model of [18] has a somewhat different linear growth rate, the structure of the <p equation is 
unchanged, and the behavior of the a = model is correspondingly very similar to that of the a — 1/2 case [18]. 
Note that while the simulated models of [18] arc similar to ours, the simulation technique is different: the simulations 
presented here solve the full <p equation (2.13), while [18] adopted an approximation neglecting the tensorial structure 
of Eq. (2.13). See Appendix A for a discussion. 

On larger scales, k ~ l/r c , the full higher-dimensional theory becomes important. Parametrized generalizations 
of DGP were proposed in [17], which lead to interesting effects on the low multipoles of the CMB. In contrast, the 
nDGP+DE model exhibits the five-dimensional features of the DGP model on large scales, except of course for the 
expansion history which is ACDM. Higher-dimensional effects however should not be important for predictions on 
sub-horizon scales. In this regime, the results presented here for the nDGP+DE model will be relevant to generalized 
braneworld models as well. 



In [1], we have introduced cosmological simulations of the self-accelerating DGP model, which self-consistently solve 
the full non-linear Eq. (2.13) of the brane-bending mode. The code, based on an earlier version presented in [30], uses 
a standard particle-mesh algorithm, augmented by a relaxation solver for <p. The Gauss-Scidel relaxation proceeds 
using Newton's method. The convergence speed is greatly enhanced using multigrid techniques [31, 32]. We use this 
code for simulations of the nDGP+DE models, adapting the expansion history and the (3 [Eq. (2.9)] function as it 
appears in Eq. (2.13). See [1] for details on the implementation and tests of the code, which all apply in an analogous 
way to the nDGP+DE simulations. 

We simulate two nDGP+DE models which only differ in their value of r c , 500 Mpc in the case of nDGP-1, and 
3000 Mpc for nDGP-2. The remaining cosmological parameters defining the expansion history and the primordial 
power spectrum are taken from the best-fit flat ACDM model given in [10]. The data used in their fit are the WMAP 
5yr results, Supernova and Ho measurements. Tab. I summarizes the cosmological parameters of the simulations. 

In addition to the full simulations, we performed simulations for nDGP-1 and nDGP-2 using the linearized field 
equation for ip, by replacing G in the GR simulations with Gu n (a) given in Eq. (2.8). We refer to the latter as linearized 
DGP simulations. We also simulated a standard ACDM cosmology with identical expansion history. Deviations of 
the nDGP+DE simulations from the ACDM simulations are thus purely due to the modification of gravity. We have 
verified that our ACDM simulations match fitting formulas for the non-linear power spectrum and halo mass function. 
Since these results are essentially the same as those shown in [33, 34], we do not show them again here. 

The simulations were started at Zi = 49. In order to test whether this is early enough to capture the non-linear 
evolution in our smallest box, we compared the power spectrum in that box after 10 time steps (a = 0.04) to the 
initial power spectrum (a = 0.02). The departures from scale-free linear evolution in the power spectrum are below 
1% even up to k = 2/c max w 6/i/Mpc. We conclude that given the limited resolution of our simulations, there is no 
need to start the simulations at an earlier time. 

However, as can be seen from Eq. (2.8), there are small residual effects from DGP even at the initial epoch, where 
(3(zi) is of order a few hundred. In order to take into account these small modified gravity effects, we correct the linear 
ACDM power spectrum at Zi as given by CAMB [35] by a factor f 2 (k) where f(k) = Tppf,dgp(&, Zi)/^ACDM(&, Zj). 
Here, 7ppf,dgp is the DGP transfer function calculated in the parametrized Post-Friedmann approximation [29], 
which has been shown to be accurate at the few percent level on large scales. On small, quasi-static scales the 
calculation is exact. Tacdm was calculated in the same way without any modification to gravity. Fig. 4 (left panel, 
solid line) shows the correction factor / obtained in this way as a function of k. At Zi — 49, the DGP effects on 
the transfer function are at the level of few 10~ 3 . Hence, a percent-level uncertainty on the correction will have a 
negligible impact on the initial power spectrum. Furthermore, for simplicity and because it is a small correction, we 
run all simulations including ACDM with the same initial conditions corrected for nDGP-1 (r c = 500 Mpc). 




(2.18) 



III. SIMULATIONS 
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k [h/Mpc] k [h/Mpc], scaled to L box = 256 Mpc/h 



FIG. 4: Left panel: Deviation of the nDGP-1 matter transfer function from ACDM at the initial redshift of the simulations, 
Zi — 49. The solid line shows the PPF calculation for DGP which is used to correct the power spectrum for the initial 
conditions. The red dashed line shows the same deviation calculated using the quasi-static calculation [Eqs. (2.8)-(2.9)]. Right 
panel: Relative deviation of the matter power spectrum at z — measured in simulations with Gaussian smoothing of the 
r.h.s of Eq. (2.13) to simulations without smoothing. The thick lines show results for Lbox = 400 Mpc/ft,, while the thin lines 
show those for Lb ox = 128 Mpc/h. All k values were rescaled to those for a 256 Mpc/h box. The vertical dashed line shows the 
maximum wavenumber used in the analysis of the simulation results, fc max = &Ny/8 (see text). 



Fig. 4 (left panel) also shows the same correction calculated in the pure quasi-static approach [Eq. (2.8)], which 
neglects the non-local Vk^/r c term in the Poisson equation [18, 22]. On scales relevant for the simulations, the 
quasi-static approximation is accurate. This is important, since the N-body simulation assumes a quasi-static regime. 
Fig. 4 shows that this is justified at Zj = 49 for the nDGP models considered here (the deviations from the quasi-static 
assumption will continue to decline at later times). 

Our simulations use a 512 3 grid with 512 3 particles. The high number of particles was chosen in order to reduce 
the shot noise in the density field. Since Eq. (2.13) is non-linear in the highest derivatives, it responds sensitively 
to small-scale inhomogeneities in the density field. Increasing the number of particles helps in reducing the residual 



TABLE I: Cosmological parameters of the simulations. 





ACDM 


nDGP-1 


nDGP-2 


— 1 f^A 


0.259 


0.259 


0.259 


r c [Mpc] 


CO 


500 


3000 







17.5 


0.487 


Ho [km/s/Mpc] 


71.6 


71.6 


71.6 


100 fi 6 h 2 




2.26 




Q c h 2 




0.110 




r 




0.0825 




n 3 




0.959 




A s (k = 0.05 Mpc -1 ) 




2.107 10" 


9 


cr 8 (ACDM) ° 


0.7892 



"Linear power spectrum normalization today of a ACDM 
model with the same primordial normalization. 



TABLE II: Simulation type and number of runs per box 
size. 





ibox [Mpc/h] 


400 256 128 64 


ACDM 


3 3 3 6 


Linearized nDGP-1 
Full nDGP-1 


3 3 3 6 
3 3 3 6 


Linearized nDGP-2 
Full nDGP-2 


3 3 3 6 
3 3 3 6 


fcmax = k Ny /8 [/l/Mpc] 


0.50 0.79 1.57 3.14 


r- C eu [Mpc/h] 


0.78 0.50 0.25 0.13 


M min [10 12 M Q /h] 


219 57.3 7.17 0.90 


r s [grid cells] a 


0.8 0.8 0.8 



"Gaussian smoothing radius for full DGP simulations. 
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errors in the numerical solution of Eq. (2.13) to an acceptable level [1]. Note that we adopt a very strict (conservative) 
convergence criterion, in demanding that the dimcnsionlcss RMS residuals of the field equation are less than 10~ 10 . 
For the self-accelerating simulations presented in [1], we employed a Gaussian smoothing of the r.h.s. of Eq. (2.13). 
Fortunately, the convergence properties of the nDGP models considered here are better, since the cross-over scale 
is smaller: this raises the density threshold for which the non-linearities in ip become important. Correspondingly, 
particle noise in the density field has less impact. This can also be seen by considering the Vainshtcin radius of a 
single particle in the simulations: for the self-accelerating model of [1], r* was of order 1 grid cell. For the models 
considered here, it is less than 0.5 grid cells. 

For this reason, we were able to eliminate the smoothing of the r.h.s. for our smallest box (64Mpc/ft. comoving 
size), where the smoothing effects are largest [1]. In addition, we lowered the smoothing radius to 0.8 grid cells for all 
other boxes. These choices were again made based on achieving an RMS residual less than 10~ 10 at all time steps. 
We quantified the effect of the smoothing by comparing full DGP simulations with and without smoothing, where the 
latter slightly violate our upper bound on the residuals (maximum residual of - 3 - 6 x 10~ 10 ). Fig. 4 (right panel) 
shows the power spectrum at z = of simulations with smoothing relative to simulations without smoothing using 
the same initial conditions. The effects of the smoothing mainly appear on wavenumbers larger than the adopted 
maximum wavenumber, given by fc max = fcN y /8, where k^ y is the Nyquist frequency of the grid [1]. Below /c max , 
the differences are below 1.5%, significantly smaller than our statistical uncertainties on the final power spectrum 
measurement (Section IV A). The same was found for the halo mass function. Hence, no correction for smoothing 
effects is necessary. 

Note that since the smoothing damps power in the brane-bending mode on small scales, and the brane-bending 
mode is attractive in the nDGP model, the simulations with smoothing of the r.h.s. have suppressed power on small 
scales. Fig. 4 also shows that the power spectrum is stable even for simulations that violate our adopted convergence 
criterion by a factor of 3-6. We simulated four different comoving box sizes, from 400Mpc/ft, to 64Mpc//i. The 
number of runs for each model and box is summarized in Tab. II. 



IV. RESULTS 
A. Matter Power Spectrum 

First, we compare the power spectrum measured in DGP simulations Pr>Gp(k) to that of the GR+DE simulations 
with the same expansion history, Pacdm(^) (Fig. 5). The procedure is the same as that employed in [1, 33]: we calcu- 
late the relative deviation of the power spectrum run by run, comparing simulations with identical initial conditions. 
We then determine the average deviation and its error by bootstrapping over realizations. In this way we are able to 
reduce the variance of the deviations considerably. 

The left panel in Fig. 5 shows the results for full and linearized simulations of the nDGP-1 cosmology at redshift 0. 
Both simulations approach the scale-invariant linear predictions on the largest scales (k < 0.1 /i/Mpc). For nDGP-1 
on non-linear scales, we find good overall agreement with the results of Khoury & Wyman (Fig. 11 in [18]). This 
is expected, since we found that the effects on the power spectrum of the G e ff-approximation used in [18] are at 
the level of 5-10 percent in our moderate resolution simulations (Appendix A). The overall scale-dependence of the 
deviations in the non-linear matter power spectrum shows the typical behavior when comparing a cosmology with 
the same linear power spectrum shape but slightly different normalization. The deviations initially grow towards 
smaller scales, peak around k ~ 0.7/i/Mpc, and then decline again towards even smaller scales. At the peak, the 
matter power spectrum in nDGP-1 is enhanced by a factor of 2 with respect to ACDM, while the enhancement is 
about 65% in the linear regime. In the halo model, the largest deviations occur at a scale roughly corresponding to 
the transition region between two-halo and one-halo terms. On smaller scales, the power spectrum mainly probes 
the density profiles of dark matter halos. If the halo profiles are not strongly affected by the increased growth, the 
deviations are expected to decrease again in this regime (see Section IV C for a study of halo profiles). These trends 
are captured qualitatively by the halof it prescription [36], which is used to map a linear power spectrum into a non- 
linear one and is calibrated on GR N-body simulations. The dashed line in Fig. 5 shows the deviation of the non-linear 
P(k) for DGP from that for ACDM obtained using halof it with the corresponding linear power spectra. While the 
halof it predictions reproduce the qualitative behavior, the overall magnitude of the deviations on non- linear scales 
is not matched. Interestingly, halof it does predict the slight suppression of the deviations below the linear prediction 
on quasi-linear scales k < 0.1 h/Mpc. A modification of the halof it parameters to improve the fit to modified gravity 
simulations was presented in [18]. Note that if there was a unique prescription to go from linear power spectrum to 
the non-linear one at a fixed redshift, then halof it should be able to describe the linearized DGP simulations, as 
the linear power spectrum in the nDGP+DE model is equivalent to a ACDM power spectrum with a slightly higher 
normalization on these scales (k > 0.01/i/Mpc). The discrepancy between the linearized DGP simulations and the 
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FIG. 5: Relative deviation in the matter power spectrum from ACDM at z = 0, AP(fc)/PACDivi(fc) = (-PdgpW — 
-PACDM(fc))/PACDM(fc) for nDGP-1 (left panel) and nDGP-2 (right panel). The blue squares and red triangles show the mea- 
surements from the full and linearized DGP simulations, respectively. The short-dashed line shows the almost scale-invariant 
deviation predicted using linear perturbation theory (Section II B), while the long-dashed line shows the deviation obtained 
when using the halof it prescription together with the linear prediction for DGP. 



halof it prediction shows that the linear-to-nonlinear mapping does depend on the growth history (which is different 
in nDGP+DE), rather than just the growth factor at the present time. 

For nDGP-1, the deviations from GR in the full simulations are only marginally suppressed compared to those in 
the linearized simulations, on the scales accessible to our simulations. This is somewhat expected, since the threshold 
overdensity for the Vainshtein mechanism in nDGP-1 is ~ 1000 (Section II C), which is only reached in the cores of 
the most massive halos in our simulations. The situation is slighly different in nDGP-2 (right panel of Fig. 5): here 
the effects of the non-linear interactions of if are noticeable for k > 0.2/i/Mpc, and become increasingly important 
towards smaller scales. The behavior as function of k of the deviation in non-linear matter power is very similar to 
that of nDGP-1, with the overall magnitude being smaller (20 — 30%). 

Note that in both nDGP cosmologies, the power spectrum deviation becomes significantly scale-dependent at around 
k ~ 0.1 h/Mpc, in the region where the second and higher harmonics of the baryon acoustic oscillations (BAO) are 
located. This scale-dependence should be taken into account when using (sufficiently precise) BAO measurements to 
constrain DGP and other braneworld models. 



B. Halo Mass Function 



Previous studies have shown that the abundance of dark matter halos is a sensitive probe of modifications of gravity 
[1, 26, 34]. We identify dark matter halos in the DGP and GR+DE simulations using a spherical overdensity halo 
finder as described in [34], with the same mass thresholds as adopted in [1], corresponding to a minimum of 6400 
particles in a halo (see Tab. II for the minimum halo mass for each box size). The halo mass M200 is defined as 
the mass contained within a spherical region of radius r2oo around the halo center of mass, whose average density is 
200 x p m . We calculate the mass function deviation of the nDGP simulations from the GR simulations run- by-run, 
and then take the average, as done for the power spectrum. 

Fig. 6 shows the relative deviation of the mass function n\ n m = dn/d\nM in DGP from ACDM, where n(M) is the 
number density of halos of a given mass. The left panel shows the result for nDGP-1, the right panel for nDGP-2. In 
each case, the top panel shows the deviation at z = 0, while the bottom panel shows z = 1. The qualitative behavior 
is the same in all cases and is expected: the increase in halo abundance grows rapidly with mass, as the abundance 
of massive halos is exponentially sensitive to the amplitude of fluctuations today, which is increased in the nDGP 
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FIG. 6: Relative deviation of the halo mass function from ACDM, An\ n m /n\ n a/(ACDM) where AninM = ni n ji/(DGP) — 
ninM(ACDM), at z — (top) and z = 1 (bottom) for nDGP-1 (7e/f panel) and nDGP-2 (right panel). Blue squares and red 
triangles denote the full and linearized DGP simulations, respectively. 

cosmologies. Since the abundance of massive halos can be probed using galaxy clusters, this observable can serve as 
a sensitive probe of braneworld gravity. Current observations should be able to put a lower limit on the cross-over 
scale r c in the Gpc range, a constraint which would be independent of the specific DGP expansion history. 

Note that at a fixed mass, the relative enhancement of the mass function at z = 1 is comparable to the one at 
z = 0. This is due to a cancelation of two effects: the growth enhancement in nDGP is smaller at z = 1 than at 
z = 0; on the other hand, a fixed mass corresponds to a rarer fluctuation at z = 1 than at z = 0, hence increasing the 
sensitivity to the growth enhancement. 

At z — 0, the Vainshtcin mechanism does not strongly affect the halo mass functions in our nDGP+DE simulations. 
For nDGP-2, there is an indication of a slightly suppressed abundance of the highest mass halos in the full simulations 
compared to the linearized simulations, which one expects since the most massive halos have the largest effective 
Vainshtein radii. This situation changes somewhat at higher redshifts: in nDGP-2 at z = 1, the mass function in 
the full simulations is lower by 5 — 10% at all masses. This can be explained by the lower threshold overdensity at 
which the (p self-coupling becomes important (Section II C): in nDGP-2, the non-linearity parameter g is ~ 0.06 at 
2 = 1, compared to 0.01 at z = 0. Since a much larger portion of the universe has S > 10 rather than S > 100, the 
self-coupling of ip affects a broad range of halo masses at z = 1 . 

C. Halo Profiles 

Once dark matter halos have been identified, we can study the distribution of mass around them. We stack the 
radially averaged density profiles of halos in mass intervals as described in [34] . The center-of-mass of each halo is 
calculated as the center of mass of particles within 1.4r co n of the center of the central (highest density) cell, which 
comprise a large fraction of the halo mass. In order to reduce scatter within the mass bin, we scale each density profile 
to its own 7"2oo before stacking. We then bootstrap over all halos in the given mass range in order to determine the 
average profile. The spatial resolution of our particle-mesh simulations is limited by the fixed size of grid cells r cc u 
(see Tab. II). We measure halo profiles down to the grid scale, though we expect that profiles have converged only at 
scales of several grid cells. We only use the highest resolution boxes (Lbox = 64Mpc//i) for the profile measurements. 

Fig. 7 (upper panels) shows the halo density profiles measured in the full DGP simulations, in three mass bins: 
10 12 - 1O 13 M //i (left), 10 13 - 10 14 M Q /h (center), 10 14 - 10 15 M Q /h (right). In the inner regions, r < r 200 , all 
DGP simulations show the same, universal profile. Moreover, the profiles match those of the corresponding GR+DE 
simulations, as shown in the lower panels in Fig. 7. This is because the inner regions of halos are assembled early 




FIG. 7: Average halo density profiles measured in the full DGP simulations. Each figure corresponds to a fixed mass range: 
lg Mj M e /h = 12- 13 (left), 13 - 14 (middle), 14 - 15 (right). The top panels show profiles of 5 P = (p- p m ) /p m vs. r in units 
of r*2oo, while the bottom panel shows the relative deviation from the GR simulation with the same expansion history in each 
case (ACDM for nDGP, QCDM for sDGP). For all figures, individual halo profiles were scaled to their respective r2oo before 
averaging. 



on in the history of the universe, where the effects of the brane-bending mode are very small (f3 3> 1). Though our 
simulations cannot probe deep into the halos, it seems unlikely there will be significant departures from GR at still 
smaller radii. 

Outside of the virial radius of halos, where matter is still infalling and tidal fields are important, modified gravity 
effects can be seen: the departures from the GR simulations typically peak around 2 — 3 r2oo , and decrease again 
towards larger radii. As expected, they are positive for the nDGP+DE simulations, largest for nDGP-1, and negative 
for sDGP. Note that the profiles for the highest mass bin are noisy due to small halo statistics (especially in the 
case of sDGP where the abundance of massive halos is significantly suppressed). For r > r 2 oo, Fig. 7 should be 
interpreted as a halo-mass correlation function (scaled to r 2 oo), which is related to the mass correlation function and 
power spectrum. Hence, the enhancement of the halo-mass correlation function in the nDGP simulations, peaking at 
a few r2oo, is another way of seeing the power spectrum enhancement, peaking at k ~ 0.7/i/Mpc (Fig. 5). 

We can also look at the behavior of the brane-bending mode within halos. Fig. 8 (top panel) shows the average 
profile of the radial gradient |Vr¥?| x 3/3/2 measured in the full DGP simulations around the most massive (and best- 
resolved) halos. We show the gradient of tp as the quantity with observable effects, since <p itself is not observable. The 
prefactor is chosen so that the quantity approaches V r <I>_ on linear scales (Section II B). Clearly, ip and its gradient 
are suppressed within halos. As expected, the suppression is most severe for sDGP, and weakest for nDGP-2. In 
each case, the thin dashed lines in Fig. 8 show the approximate solution for ip presented in [18], calculated given the 
density field as described in [1] (see also Appendix A). While the approximation works well at the very center of the 
halos, it does not describe the solution at r > 0.5 r 2 oo- Moreover, the discrepancy becomes larger the stronger the 
non-linearities in the ip equation are, i.e. the larger r c is. 

The lower panel of Fig. 8 shows the gradient of the dynamical potential around the same halos compared to the 
gradient of the Newtonian or lensing potential = which obeys the usual, unmodified Poisson equation. A 
discrepancy between the dynamical potential of a given mass and its lensing potential is precisely the quantity that 
is constrained by Solar System tests, described via the post-Newtonian parameter 7 [37]. Fig. 8 shows that this 
deviation is indeed suppressed within massive halos in our simulations, although the simulations lack the resolution 
to follow the suppression very far into the cores of halos. As expected, the suppression is strongest for sDGP which 
has the largest non-linearity parameter g (Section II C), and weakest for nDGP-1. 



D. Bispectrum 



The matter bispectrum in DGP is, in principle, a probe of the quadratic non-linearity in the brane-bending mode 
equation. This non-linearity adds an additional contribution to the tree-level bispectrum generated by the ordinary 
gravitational non-linearities [22]. As discussed in Section II C, the specific quadratic self- interactions of the brane- 
bending mode lead to a dependence on geometry. In Fourier space, the non-linearity in Eq. (2.13) leads to a coupling 
between two Fourier modes <£(ki), <^(k 2 ) with a kernel proportional to 1 — (k\ ■ k^) 1 , where ki denote unit vectors of 
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FIG. 8: Top panel: Average profiles of the radial gradient of ip, i.e. the force modification, measured in the full DGP simulations 
and scaled by 3/3/2, for halos of mass lg Mj Mq/H = 14 — 15. The profiles were averaged in a similar way as the density profiles. 
Bottom panel: Deviation of the dynamical potential ^ [Eq. (2.10)] from the Newtonian potential or lensing potential $_ = ^jv 
[Eq. (2.12)] around the same dark matter halos. This quantity is probed by Solar System tests. 




FIG. 9: Ratio of the reduced bispectrum measured in full and linearized DGP simulations, vs the angle between the wavevectors 
ki = 0.5 fa/Mpc, k2 = 1 /i/Mpc. This is a direct probe of the configuration dependence of the Vainshtein effect. For each model, 
the bispectra were measured in the 256Mpc/fo boxes. The errors are bootstrap error bars on the ratio determined from the 
different runs. We have omitted the error bars on the nDGP-2 model for clarity - they are very similar to those for nDGP-1. 



the respective wavenumbers. As shown in [22], this structure manifests itself as a characteristic change of the shape 
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dependence of the bispectrum, which is commonly expressed in terms of the reduced bispectrum Q: 

ni\r v n = g(ki,k 2 ,k 3 ) 

Q[ U 2 ' 3) " P(kr)P(k 2 ) + P(k 2 )P(k 3 ) + P(k 3 )P(k 1 ) ' (4 ' ij 

where P(ki,k 2 ,k 3 ) is the bispectrum. In order to separate these characteristic effects of the brane-bending mode 
interactions as cleanly as possible from the ordinary gravitational non-linearities, we compare the bispectrum measured 
in the full DGP simulations to that of the linearized DGP simulations, which have an identical linear growth factor, 
and hence are very similar in their state of non-linear evolution. Fig. 9 shows the ratio of the reduced bispectrum 
measured in the two cases, for fixed wavenumbers of k\ = 0.25 /i/Mpc and fc 2 = 0.5 /i/Mpc, as a function of the angle 
between them, cos = k\ ■ fc 2 • The bispectrum was measured in the simulations using Monte Carlo integration in a very 
similar way as described in [38]. We again average the ratios run-by-run, and estimate error bars by bootstrapping 
over realizations. 

The strongest effects are seen for the sDGP case, since this cosmology has the largest r c and hence the lowest 
threshold for the onset of the ip self-interactions. Note that the bispectrum is enhanced in the full sDGP simulations 
for equilateral configurations 9 « 7r/2, while it is very close to the linearized simulations for squeezed configurations 
6 = 0, it. For squeezed configurations, which correspond to planar geometry, the kernel 1 — {k\ ■ fc 2 ) 2 = and the 
non-linearities vanish. For equilateral configurations, the suppression of ip is strongest, which leads to an enhanced 
bispectrum in the full DGP simulations since the brane-bending mode is repulsive. For nDGP-1 and nDGP-2, we 
see a hint of the opposite effect: here <p is attractive, and a suppression due to the self-interactions in equilateral 
configurations leads to a suppression of the bispectrum relative to the linearized simulations. The overall effect is 
much smaller, since r c is smaller and the non-linearities in ip only become important at higher densities (i.e., smaller 
scales). 

The result for the sDGP bispectrum agrees very well with that found in [26]. Note that lacking very large box sizes, 
we measure the bispectrum at relatively small scales. The ordinary gravitational non-linearities which enter at higher 
order already contribute at the percent level on the scales accessible in our simulations. Hence, we do not attempt a 
quantitative comparison with tree-level perturbation theory predictions here. 



V. DISCUSSION 



We have introduced the nDGP+DE model based on the normal branch of DGP with a smooth dark energy 
component on the branc, which results in an expansion history that is precisely ACDM. This model can serve as an 
approximation, especially on sub-horizon scales, to more general braneworld models whose cosmological solutions have 
not been obtained yet. In addition, the r c — > oo limit precisely corresponds to General Relativity plus cosmological 
constant. 

Geometric measurements such as the acoustic scale constrained from the CMB, BAO, and Supcrnovae observations 
do not constrain the crossover scale r c between 4D and 5D gravity in this scenario. The growth of large-scale structure 
however does offer a sensitive probe of braneworld gravity via the effects of the brane-bending mode which mediates 
an additional attractive force. In the case of the self-accelerating branch, the brane-bending mode is repulsive, so that 
the modified gravity effects are sign-flipped. 

In order to study the formation of structure in these models, it is necessary to solve the non-linear equations of 
the brane-bending mode in conjunction with the ordinary gravitational dynamics in N-body simulations [1]. These 
non-linear interactions, which are responsible for the Vainshtein mechanism [23, 24] restoring General Relativity in 
high density environments, also leave characteristic signatures in the bispectrum in the simulations [22], which we 
have found in both sDGP and nDGP simulations. In addition, the non-linearities manifest themselves in a suppression 
of the brane-bending mode around massive halos. 

The matter power spectrum in nDGP+DE shows a characteristic enhancement compared to that of ACDM, in- 
creasing towards smaller scales up to k <~ 0.7ft./Mpc, and decreasing on even smaller scales. This is evidenced in 
a complementary way by the profiles of dark matter halos and their environments (halo-mass correlation function). 
Similar to what was found in f{R) simulations [34], the inner regions of halos are not affected by the enhanced forces 
since they assembled at early times. The strongest effects are seen in the near environment of halos, around 2 — 3 r 2 oo, 
corresponding to the scale where the power spectrum enhancement peaks. 

The abundance of massive halos is known to be a sensitive probe of the growth of structure [34, 39-41]. We found 
that indeed massive clusters are 2 — 4 times more numerous for the nDGP+DE model with r c = 500 — 3000 Mpc. 
Current cluster abundance measurements (e.g. [42-44]) should be able to constrain the crossover scale r c to be at 
least of order Gpc - a constraint which can be placed independently of the specific DGP expansion history. 

As the next step in understanding the effects of braneworld gravity on large scale structure, an upcoming paper 
will detail a model of the matter power spectrum and halo mass function based on spherical collapse and the halo 
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model. Such a model can also serve as an efficient way of dealing with cosmological parameter dependencies when 
comparing with actual observations [27]. In the future, by means of these simulations and a physical model of their 
results, we will hopefully be able to probe the next generation of braneworld scenarios via their effect on the growth 
of cosmic structure. 
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APPENDIX A: COMPARISON WITH G efT (S) APPROXIMATION 

We now discuss the approximation adopted in the simulations of [18] and compare it with the results of our full 
solution of the brane-bending mode equation (2.13). Consider a spherical "tophat" mass with a constant density 
contrast 5p. Then, the ansatz ip(r) = Ar 2 + C solves Eq. (2.13), and in particular we have: 

(V«V^) 2 = i(vV) 2 - (Al) 

For this special case, the non-linear terms combine to an equation of ip which is algebraic in V 2 <p: 

„■> 2r r n2 8irGa 2 . , . „. 

vV+ 9^ (vV) =-3r*»* (A2) 

This can then be solved for V 2 ip in terms of a non-analytic function of the overdensity 5: 

VV = 8nG cB (S)a 2 p m S, (A3) 



G cff = ^ , 3 5. (A4) 

Thus, the brane-bending mode is determined in this approximation by a Poisson equation with density-dependent 
gravitational constant G e s(S). Since this is a linear equation in ip which can be solved via Fourier transform, a 
simulation using this approximation is similar in terms of computing time as ordinary GR simulations, and much less 
computationally expensive than solving the full non-linear differential <p equation. 

We now discuss the caveats of this "G e ff approximation". For more general density profiles 5p ^ const., Eq. (A3) is 
only an approximate solution. More importantly, as already pointed out in [18], the approximation does not have the 
correct large-distance behavior even for an isolated top-hat mass, so that scales in the linear regime of cosmological 
perturbations can in principle be affected by this incorrect behavior. 

To see this, consider an isolated mass with dp = for r larger than some radius R. Integrating Eq. (2.13) over a 
sphere with radius r > R, we see that the two non-linear terms in Eq. (2.13) cancel via partial integration, leaving only 
boundary terms which become increasingly suppressed as we let r — > oo. Since the right-hand side is just proportional 
to the enclosed mass M within r, we see that ip approaches the linearized solution ip(r) = (2/3(3)GM/r at large r, 
irrespective of the strong non-linearities close to the mass. For this reason, we recover the linear predictions of the 
DGP model (Section II B) on large scales in the simulations. 

In contrast, integrating Eq. (A3) over the sphere with radius r > R, we see that the asymptotic behavior in this 
approximation is ip(r) — > (2/3f3)G e gM/r, where G c g is the effective gravitational constant averaged over the mass. 
As Eq. (A4) shows, G c g <C G if the density contrast 5 ^> 1. Hence, the far field of <p is strongly suppressed in this 
approximation. In particular, if most of the cosmological mass is in non- linear structures with S > 1, we do not expect 
to recover linear theory predictions on large scales in this approximation. 

As this discussion shows, the error encurred with the G e s(6) approximation in cosmological N-body simulations 
depends on the resolution of the simulations: very low resolution simulations (with large box sizes) which resolve 
little non-linear structure will not be affected; high-resolution simulations which resolve very small dark matter 
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FIG. 10: Relative deviation of the power spectrum at z = measured in a simulation using the G e g(5) approximation, PG a!f , 
from that of the full DGP simulation, PdgpW f° r nDGP-1 and nDGP-2. In each case, results are shown for a "typical" run 
with L box = 64Mpc/fr. 



halos will show a strong suppression of the <p field and the corresponding modified forces on large scales, if the 
G e fi approximation is used. In order to test how big the effect is in our moderate resolution simulations, we reran a 
representative simulation of our smallest box (Lbox = 64Mpc//i) for the nDGP-1 and nDGP-2 models using Eq. (A3) 
instead of Eq. (2.13). For this, we chose a realization with a power spectrum close to the average of our 6 realizations 
of the 64Mpc//i box. Fig. 10 shows the relative deviation in the power spectrum at z = in simulations with the 
G e ff approximation from that of the full DGP simulation with the same initial conditions, for nDGP-1 and nDGP-2. 
The deviations are around 5-10% for nDGP-1 and ~ 3% for nDGP-2. The magnitude of deviations is noticeable 
given our < 1% precision on the power spectrum deviation on quasilinear scales (Fig. 5). As expected from our 
discussion, the deviations persist on the largest scales probed by this box size. Furthermore, power is suppressed in 
the G e ff simulations due to the suppression of the attractive (/^-mediated force. The deviations are larger in nDGP-1 
simply due to the stronger effect of ip on structure formation in the small-r c model. Thus, while the use of the G c s(8) 
approximation would not have a dramatic impact on the main results presented in this paper, care must be taken in 
future high-resolution simulations of braneworld models: the artefacts of this approximation will increase rather than 
shrink with increasing resolution. 
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